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Abstract 

Poland-Scheraga models were introduced to describe the DNA denaturation transition. We 
give a rigorous and refined discussion of a family of these models. We derive possible scaling 
functions in the neighborhood of the phase transition point and review common examples. 
We introduce a self-avoiding Poland-Scheraga model displaying a first order phase transition 
in two and three dimensions. We also discuss exactly solvable directed examples. This 
complements recent suggestions as to how the Poland-Scheraga class might be extended in 
order to display a first order transition, which is observed experimentally. 

1 Introduction 

When a solution of DNA is heated, the double stranded molecules denature into single strands. In 
this process, looping out of AT rich regions of the DNA segments first occurs, followed eventually 
by separation of the two strands as the paired segments denature. This denaturation process 
corresponds to a phase transition |35|. 

A simple model of the DNA denaturation transition was introduced in 1966 by Poland and 
Scheraga '30' "ST (hereinafter referred to as PS) and refined by Fisher TD' TT. The model consists 
of an alternating sequence (chain) of straight paths and loops, which idealize denaturing DNA, 
consisting of a sequence of double stranded and single stranded molecules. An attractive energy 
is associated with paths. Interactions between different parts of a chain and, more generally, all 
details regarding real DNA such as chemical composition, stiffness or torsion, are ignored. It was 
found that the phase transition is determined by the critical exponent c of the underlying loop 
class. Due to the tractability of the problem of random loops, that version of the problem was 
initially studied by PS 31 . The model displays a continuous phase transition in both two and 
three dimensions. It was argued ,10. that replacing random loops by self-avoiding loops, suggested 
as a more realistic representation accounting for excluded volume effects within each loop, sharpens 
the transition, but does not change its order. 

However, the sharp jumps observed in the UV absorption rate in DNA melting^ experiments 
|35| . which correspond to a sudden breaking of large numbers of base pairs, indicated that a 

^In recent years, a number of other properties of DNA molecules have been studied by refined techniques such as 
optical tweezers and atomic force microscopy, and theoretical descriptions of underlying efli'ects such as unzipping 
1211 1:^31 1^1^71 have been proposed. These will not be discussed here. 
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first order phase transition would be the appropriate description. The question whether such an 
asymptotic description, which imphes very long chains, is valid for relatively short DNA sequences, 
has been discussed recently 22 . Nevertheless, a directed extension of the PS model, being 
essentially a one-dimensional Ising model with statistical weighting factors for internal loops, is 
widely used today and yields good coincidence of simulated melting curves with experimental 
curves for known DNA sequences |85l HI IK . Another recent application of PS models analyses the 
role of mismatches in DNA denaturation 134. A numerical approach to DNA denaturation, which 
we will not discuss further, uses variants of the Peyrard-Bishop model |29l I84j . a Hamiltonian 
model of two harmonic chains coupled by a Morse potential. 

With the advent of efficient computers, it has more recently been possible to simulate analyt- 
ically intractable models extending the PS class, which are assumed to be more realistic repre- 
sentations of the biological problem. One of these is a model of two self-avoiding and mutually 
avoiding walks, with an attractive interaction between different walks at corresponding positions 
in each walk The model exhibits a first order phase transition in d = 2 and d = 3. 

The critical properties of the model are described by an exponent c' related to the loop length 
distribution |20[ I21L 151 151 fT|. see also Fisher's review article ^T]. This exponent is called c 
again. Indeed, for PS models, it coincides with the loop class exponent c if 1 < c < 2, see below. 
Within a refined model, where different binding energies for base pairs and stiffness are taken 
into account, the exponent c' seems to be largely independent of the specific DNA sequence and 
of the stiffness of paired walk segments corresponding to double stranded DNA parts (S] . There 
are, however, no simulations of melting curves for known DNA sequences which are compared to 
experimental curves for this model. 

An approximate analytic derivation of the exponent related to this new model was given by 
Kafri et al. [2011^123 using the theory of polymer networks. They estimated the excluded volume 
effect arising from the interaction between a single loop and two attached walks. This approach 
(refined recently '3') yields an approximation to the loop length distribution exponent c', which 
agrees well with simulation results of interacting self- avoiding walk pairs |HlEm- There is a recent 
debate about the relevance of this approximation to real DNA ^| 122] • The polymer network 
approach, as initiated previously PDI, led to a number of related applications |21 |2S1 El 121 E| ■ 

In this article, we reconsider PS models for three principal reasons. Firstly, the older articles 
are short in motivating the use of particular loop classes, which may be misleading in drawing 
conclusions about the thermodynamic effects of different loop classes. In fact, the loop classes 
discussed in the early approaches [M\ ITHj are classes of rooted loops and lead to chains which are 
not self- avoiding. This seems unsatisfactory from a biological point of view, since real DNA is 
self-avoiding. Secondly, the common view holds that PS models with self-avoiding loops cannot 
display a first order transition in two or three dimensions. In fact, this view led to extending the 
PS class 1^1 1201 121L 1231 in order to find a model with a first order transition. However this view 
is incorrect, as we demonstrate by a self-avoiding PS model with self-avoiding loops. Thirdly, the 
two exponents c and c', extracted from different expressions as described above, are used in the 
literature without distinction, although there are subtle differences, which we will point out. 

This paper is organized as follows. In the next section, we give a rigorous discussion of PS 
models and their phase transitions. We derive the scaling functions which describe the behav- 
ior below the critical temperature, as the critical point is approached. We then derive the loop 
statistics, thereby analyzing the occurrence of the loop class exponent in the loop length distri- 
bution. The third section reviews the prevailing PS models, with emphasis on motivation for 
the underlying loop classes. We then introduce a self-avoiding PS model displaying a first order 
phase transition in two and three dimensions. In Section |31 exactly solvable directed PS models 
are discussed. We will give explicit expressions for generating functions. The critical behavior of 
some of these models has been analyzed previously by different methods |2S1 12()l ITfj . This is 
followed by a discussion of models extending the PS class. We conclude with a discussion of some 
open questions. 
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2 Poland- Scheraga models: general formalism 



We define PS models and consider analytical properties of their free energies, thereby refining 
previous expositions (301 1^ ^1 in a rigorous manner. We analyze phase transitions, extract 
the asymptotic behavior of chains at temperatures near the phase transition by means of scaling 
functions, and consider the loop length distribution of chains. 

2.1 Definition of the model 

We consider a discrete model, defined on the hypercubic lattice TU^ . Double stranded DNA seg- 
ments are modeled by paths, and single stranded DNA segments are modeled by loops on the 
edges of the lattice. Each loop is assumed to have two marked vertices to indicate where paths 
are attached. (Like PS, we do not a 'priori assume the "DNA condition" that the marked ver- 
tices divide the loop into parts of equal lengths. There is a recent discussion about the effect of 
mismatches |13|.l Any alternating sequence of paths and (marked) loops is called a chain. A PS 
model consists of all chains obtained by concatenation of paths and loops from a given path class 
and a given loop class, where the initial segment and the final segment of a chain are both paths. 
(The case of an open end affects the behavior above the critical temperature and will not be 
considered here Note that, in general, such chains are not self- avoiding, in contrast to real 
DNA. Self-avoidance may be violated by paths, by loops, or by the way segments are concate- 
nated. A chain is called segment-avoiding if there are no overlaps, i.e., every two non-neighboring 
segments have no vertex in common, and every two neighboring segments have exactly the marked 
vertex in common. We call a PS model self-avoiding if paths and loops are self-avoiding and if all 
chains of the model are segment-avoiding. We will discuss below several examples of self-avoiding 
PS models, which arise from self-avoiding walks. 

The requirement of self- avoidance restricts the admissible path classes and loop classes. A 
simple subclass of self-avoiding PS models are directed PS models: We call a chain directed if 
there is a preferred direction such that the order of the chain segments induces the same order 
on the vertex coordinates (w.r.t. the preferred direction), for each pair of vertices taken from two 
different chain segments. Such chains arc then segment-avoiding. We call a PS model directed 
if paths, loops and chains are directed. We will discuss below several examples of directed PS 
models which arise from classes of directed walks. 

For a given PS model, let z^.n denote the number of chain configurations with m contacts and 
length n. The generating function is defined by 

Z{x,w)=Y,Zm,nW"'x^, (2.1) 
rn.n 

where the activity x is conjugate to the chain length n. The Boltzmann factor w = e~^/^^ takes 
into account the attractive interaction (achieved by setting the energy £' < 0) between bonds. T is 
the temperature, and k is Boltzmann's constant. Over the relevant temperature range < T < oo, 
we have oo > w > 1. Note that there is no interaction between different segments in a chain. 

The generating function Z{x,w) can be expressed in terms of the generating functions for 
paths V{x) and loops U{x). These are 

oc oo 

V{x) = J2 ^(^) = J2 P2na;", (2.2) 

n— n— 1 

where &„ is the number of paths of length n, and p2n is the number of loops^ of length 2n. Due 
to the chain structure, we get a geometric series in V{wx)U{x), 

'M-^^^^^^-tz,.Hx-. (2.3) 

\ / \ / n— 1 

^In simulations of melting curves like MELTSIM |^ a cooperative parameter a ^ 10~^ is used. In our 
context, this amounts to replacing the loop generating function U(x) by aU{x). 
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Since we want to analyze phase transitions of the model, which can only occur in the infinite 
system, we define the free energy of the model as 

f{w) = lim - log Zn{w) = - log XcH, (2.4) 

n^oc n 

where, for fixed w, Xc{w) is the radius of convergence of Z{x,'w). Concatenation arguments and 
supermultiplicative inequalities can be used to show that the free energy exists JHl- We will 
alternatively investigate properties of the free energy in terms of the generating functions for 
paths and loops'^. Throughout the paper, we employ the following assumption. 

Assumption 1. Assume that the generating functions U{x) and V{x), defined in ^2.^) . have 
radius of convergence xjj and xy , respectively, where < xjj < xy < 1- At the critical point 
X — Xy , assume that V{xy) := lim^^^- V{x) = oo. 

Remark. U{x) is a generating function (i.e., a series with non- negative coefiicients) , such that 
U{x) and the derivative U'{x) are strictly positive for < a; < xu- A corresponding statement 
holds for V{x). The assumption that < xu < 1 reflects the requirement that the number of 
configurations p2n grows exponentially in length. The models commonly discussed PHI IS El d 
nT\ have Xy — 1 and V{0) = 1. For loops and paths of the same type, we have xjj = Xy < xy < 1, 
see also the examples discussed below. The assumption V{xy) = oo is satisfied for typical classes 
of paths. Note that typically U{x^) := lim^_^^- U{x) < oo, but there are loop classes where 

U{x~[j) — oo such as convex polygons 0. We have U{Q) =0 by definition. 

The radius of convergence Xc{w) of Z{x,w) is the minimum of the radius of convergence of 
U{x)V{wx) and the point xi{w) where the denominator in l|2.3|l vanishes. Define F{x,w) := 
U {x)V {wx) . As a function of argument x, F{x, w) is continuous and monotonically increasing for 
< X < miii{x u,xv/w). Note that F{0,w) = 0. If w > xy/xu, then VluIx^xv /w Fix^w) — oo, 
such that there exists a unique solution xi{w) < xjj with F{xi{w),w) — 1. Assume that 1 < 
w < xy/xjj. If U{x^) = oo, there exists a unique solution Xi{w) < xij. If U{x^) < oo and 
U {x'[j)V {xjj) > 1, we have F(xij,w) > 1 for all w > 1, i.e., there exists a unique solution xi{w) < 
Xjj. If U{x^) < OO and U{x^)V{xu) < 1, we define Wc > 1 by the condition U{xj^)V(wi xjj) = 1. 
Then, for w > Wc, we have a unique solution xi{w) < xu- If w < Wc, there is no such solution, 
and the dominant singularity of Z(x, w) occurs at a; = xjj. We thus proved the following theorem. 

Theorem 1. If Assumption 1 is satisfied, the PS model has free energy 

f(u,)^l -log:^c(^) {w,<w <oo) 
^ ^ \ -logxu {Kw < Wc), ^ ' 

where the radius of convergence Xc{w) of Z{x, w) is, for Wc < w < oo, the unique positive solution 
of 

U{xc{w))V{wXc{w))^l. (2.6) 

IfU{xjj) := lini^^^- U{x) = oo or U{xjj)V{xu) > 1, we have Wc — 1. Otherwise, Wc > ^ is 
implicitly given by U{x^)V{wcXu) ~ I. □ 

Remark. If Wc > 1, this point is a critical point, see the following section. 
2.2 Phase transitions 

For w > Wc, there is no phase transition, since the radius of convergence Xc{w) > is analytic 
for < w < oo, as can be inferred from H2.6|l . A necessary condition for a phase transition at 
a finite temperature > 1 is U{x^)V{xij) < 1. The nature of the transition at ui = Wc > 1 is 

^The reader may find it illuminating in following this very general discussion that now follows to refer to a 
specific model, discussed in Section Fs.ll 
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exhibited by the fraction of shared bonds 6{w) = w j^' , which is defined for w ^ Wc- Note that 
9{'w) — for w < Wc- For w > Wc, we have 

Consider the hmit w We have < V < oo, since UV = 1 and U is finite. Since V{x) is a 

generating function, this impUes < V < oo. If U'{x^) = oo, it follows that 0{w) 0, such that 
the phase transition is continuous. If U'{x^) < oo, 9(w) — > 6'c > 0, such that the phase transition 
is of first order. This leads to the following statement. 

Theorem 2. Let Assumption 1 he satisfied. IfU{x^) = oo or U{x^)V{xu) > 1, the PS model 
\2.!^) has no phase transition at finite temperature. Otherwise, ifU{xj}) < oo and U'{xj}) — oo, 
a continuous phase transition will occur at w ~ Wc > ^ defined in Theorem 1. IfU{x^) < oo and 
U'{x^) < oo, a first order phase transition will occur at w = Wc > 1 defined in Theorem 1. □ 

Remark. The phase transition condition U{x'[j)V{xu) < 1 is typically satisfied in more realistic 
models, where U{x) is multiplied by the cooperativity parameter a ~ 10~^. 

We conclude that the nature of the transition is determined by the singularity of the loop 
generating function U{x) at a; = xij. It can be more directly related to the asymptotic properties 
of loops, which are typically of the form^ 

P2n Ax'^/'n-'' {n oo), (2.8) 

for some constants A> Q and c e R. The exponent c determines the singularity a,t x — xjj, which 
is, to leading order and for c e (Q) \ N, algebraic. This specializes Theorem 2. 

Proposition 1. Let Assumption 1 he satisfied and assume that, for c e M \ N, the singular part 
of U{x) is, as x ~f xjj , asymptotically given by 

U^''"a\x) ^ Uoixu - xy-^ {x^x^) (2.9) 

for some constant Uq ^ 0. Then, the PS model 1^2. 3\) has no phase transition if < c < 1 or 
if U {x^)V {xjj) > 1. Otherwise, a continuous phase transition occurs for 1 < c < 2, and a first 
order transition occurs if c> 2, at w = Wc > i defined in Theorem 1. □ 

Remark. If c = 1 in (|2.8I) . we get a logarithmic singularity in U{x), hence no phase transition. If 
c = 2 in H2.8|l . we get a logarithmic singularity in the derivative U'{x), resulting in a continuous 
phase transition. See also Fisher 

2.3 Scaling functions 

In the vicinity of a phase transition point, critical behavior of the form 

Z{x,w) {xc - x)~^ F{{w - Wc)/{xc - x^) {x,w) ^ {Xc ,Wc), (2.10) 

uniformly in w, is expected with critical exponents 9 and (j). Here F(s) is a scaling function which 
only depends on the combined argument s = {w — Wc) / {xc — x)'^ . The scaling function is extracted 
from the generating function Z(x, w) by replacing w using the variable s of combined argument 
and expanding to leading order in x^ — x. To this end, let us assume as in H2.9(l that U{x) behaves 
as 

U{x)^U{x-) + Uo{xu -xY~^ +o{{xu -xf-^) [x^x^), (2.11) 

and 1 < c e M \ N. In this case, the constant Uq is negative. (If c < 1, no phase transition occurs 
at positive temperature according to the considerations above). We have the following result. 

^A„ ^ Bn for n ^ oo means that lim,j^oo A„/Bn = 1. Similarly, f(x) ~ g{x) for x ^ Xc means that 
limx^xe f{x)/g{x) = 1. 
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Theorem 3. Let Assumption 1 he satisfied. If the leading singularity of the loop generating 
function U{x) is of the form ^2.11}) with an exponent c > 1 and if U {x~^)V {xjj) < 1, the PS model 
h2. has critical exponents and scaling function f2.1U\l 



f, = c-l, F{s) = — (l<c<2), 

\Uo\ - U"^{x^)V'{wcXu)xus 

= 9^=1, F{s) = :: (c> 2), 

U'^{Xjj)V'{WcXu)[Wc - xus] 



(2.12) 



where Wc > 1 is implicitly given by U{x^)V{wc xu) = 1. □ 

Remark. If c > 2, the critical exponents are independent of c. In both cases, the scahng function 
is proportional to Fa{s) = 1/(1 — as) with a positive constant a. 

It can be shown jl8| that under mild assumptions such a scaling function implies a certain 
asymptotic behavior of the function Zn{w) for large n. That is to say, 



(2.13) 



uniformly in w, where h{x) is the finite-size scaling function h{x) — X^fc^o fkx'' /^{k(f> + 9), where 
T{z) denotes the Gamma function, and the coefficients fk appear in the Taylor expansion of the 
scaling function F{s) — X]fc!=o fk^''- This can be used to derive the critical behavior of the fraction 
of shared bonds above the phase transition ^H] for 1 < c < 2, 

0{w) ^ A{w - Wc)^ (u; ->«;+), (2.14) 

for some constant A > 0. In our case Fa{s) = 1/(1 — as), we have ha{x) — e'^^, ii 9 — 1. For 
O<0 = c— 1<1 and 9 rational, the finite size scaling function can be expressed in terms of 
hypergeometric functions. For the case c = 3/2 = 9+1, which will be relevant below, we get 

ha(x) = ^ + axe(°^)' (1 + erf(a2:)) , (2.15) 

where erf(a;) = 2/y^ e^* dt denotes the error function. 

2.4 Loop statistics 

Within a given PS model, consider the set of all chains of length n with m contacts. Denote the 
number of loops of length 21 in this set by gm,n,i- Define the generating function 

L{x, w,y)=Y. 9m,n,i w'^x^yK (2.16) 

Again, due to the directedness of the model, this generating function can be expressed in terms of 
the loop generating function JJ = U{x), U = U{xy) and the path generating function V = V{xw). 
An argument as in Section [TTl vields 

L{x, w, y) = VUV + {VUVUV + VUVUV) + ■■■ 

= VUV{1 + 2UV + 3{UV)^ + ■■■) ^2 17) 

= ^^^JT^JJVf ^ U{yx)Z\x,w). 

The generating function for the total number of loops is then given by setting y = 1, 

T{x, w) = L{x, w, 1) = U{x)Z'^{x, w), (2.18) 
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and the generating function for the sum of all loop lengths is given by 

= xU'{x)Z^{x,w). (2.19) 



S{x,w) = y-^L{x,w,y) 



The finite size behavior of these quantities about the phase transition can be computed from the 
scaling function results of Section [TSl For the total number of loops, we get 

\im [x"]T{x,w) ^ Ax^"-^^n^'^-^ {n~^oo), (2.20) 

where A > is some amplitude^, and for the number of loops of length 21^ we get 
lim [x"y']L{x,w,y) ^ pW^i^-^'n^'-^ (n ^ oo) 

« A r"a;^"-2^n2''-i (1< ? < n, n ^ c5o). 



(2.21) 



Note that we obtained (|2.21(l under the assumption that I is asymptotically large but that I <^n. 
This is due to the estimate [y']C/(ya;) = p2i x^ ^ P2ix\j for x — > a;^ and / finite, but p2i ^ Axj^H^'^ 
for Z ^ oo, where we assumed (|2.8|) to be valid. 

For the sum of the loop lengths, we distinguish between the cases 1 < c < 2 and c > 2, where 
the derivative of the loop generating function U'{x^) is finite or infinite at the critical point. We 
find 

hrn[x-]S{x,w) ^ \ ^1^-^^' ^ 9^ i^^^)- (2-22) 

Let us compute the probability of a loop of length 21 within chains of length n. We get, in the 
limit w — > w'^ , for the loop length distribution 

P{l,n)= lim i^y'\^(^'^^y) ^Al-'^- (i«?«n, n^oo). (2.23) 
m^wt [x'^\T{x,w) 

Furthermore, for the mean loop length ( I )„ we get, in the limit w —* w^, 

[x'']S{x,w) r An^-'^ (1< c< 2) 



(=0 



^ w^w+ [x'^\T[x,w) [An (c> 2) 



The behavior of the mean loop length refiects the nature of the phase transition: If 1 < c < 2, the 
mean loop length diverges, such that the transition is continuous. If c > 2, the mean loop length 
stays finite, indicating a first order transition. 



3 Two prominent loop classes 

Since the critical behavior of PS models is essentially determined by the properties of loops, PS 
inni 1^ , and later Fisher JHI , were led to consider various loop classes (together with straight 
paths for the double stranded segments). Whereas PS analyzed loop classes derived from random 
walks. Fisher considered loop classes derived from self-avoiding walks. 



3.1 Loops and walks 

An oriented, rooted loop of length n is a walk of length n — 1, whose starting point and end point 
are lattice nearest neighbors. As usual |28[ Sec. 3.2], we identify such loops if they have the same 
shape, i.e., if they are equal up to a translation, possibly followed by a change of orientation. These 
objects we call unrooted, unoriented loops, or simply loops. Each loop of length n has at most 2n 

^In order to simplify notation, we will denote all following amplitudes by the letter A, with the convention that 
their values may be different. 
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corresponding walks. If the walks are self-avoiding, each loop has exactly 2n corresponding walks. 
The number of loops of length n is denoted by p„. For example, for self-avoiding loops on we 
have p4 = 1 and pg = 2. For a given class of walks, the above description defines the corresponding 
(unmarked) loop class. Within a chain structure, two paths are attached to each loop. Different 
choices for attachment positions increase the number of (marked) loop configurations to pn > Pn ■ 
If we assume the DNA condition that the two paths attached to a loop bisect it into pieces of 
equal length, then the number of possible attachments of two paths to a loop of length n is less 
than or equal to 2n. (We distinguish the two strands of a marked loop). PS and Fisher consider 
classes of oriented rooted loops. Self-avoiding oriented rooted loops can be interpreted as loops 
with 2n possible attachment positions of paths to a loop of length n. A similar interpretation for 
oriented rooted random loops is not obvious. We stress that these loop classes result in chains 
which are not segment-avoiding, as paths will intersect the loops. Both models cannot therefore 
represent real (self-avoiding) DNA. 

3.2 Oriented rooted random loops 

The first simple example of loops, which are discussed by PS are oriented rooted random 
loops derived from random walks. Random walks on 'E'^ have the generating function Vd{x) = 
1/(1 — 2dx). The asymptotic behavior of the number of oriented rooted random loops of length 
2n is given [23, App. A] by 

P2„~A(2d)2"(2n)-^/2 {n^oo). (3.1) 

This implies no phase transition in c? = 2 and a continuous phase transition in d = 3, since the 
phase transition condition U{x'^)V{xu) < 1 is satisfied in d = 3: We have 

P^-- E (fc.p mpL.p - 3f2(-n, -n, -n; 1, + 1/2; 1/4), (3.2) 

k+l+m=n ^ V \ V V V y 

where 3^2(01, 02, 03; 6i, 62; 2) is a hypergeometric function. With xu = 1/36, we extracted 
the amplitude U{x^) numerically '14' using the values p2n where n < 30, getting U{x^) — 
0.51638461326(7). The result follows for straight paths since 1/(1 — xu) — 36/35 and also for 
random walks since V^{xu) — 6/5. 

3.3 Oriented rooted self-avoiding loops 

The PS results led to the question ^HI whether accounting for excluded volume effects within a 
loop increases the loop class exponent c, which might change the order of the phase transition. 
This led to considering self-avoiding loops, which are loops derived from self-avoiding walks. (A 
self-avoiding walk Sec. 1.1] is a random walk which never visits a vertex twice.) By definition, 
the self-avoiding loop class fully accounts for excluded volume interactions within a loop. In d = 3, 
self-avoiding loops of length n > 24 may be knotted. Self-avoiding walks and loops are well studied 
objects [23 ^l- Fisher considered oriented rooted loops p2n = ^np^n- Their loop class exponent 
c, P2n ~ Bii^n~'^, is related to the mean square displacement exponent v of self-avoiding walks 
by the hyperscaling relation |28[ Sec. 2.1] c = dv, where = 1/2 for d > 4, = 0.5877(6) for 
d = 3 and v = 3/4 for d = 2. At present, there is no proof for the values of v in dimensions d < 4. 
Explicitly, we have 

( d/2 d>4 
c= < 1.7631(18) d = 3 (3.3) 
[3/2 d=2 

For unknotted self-avoiding loops, which is the preferrable model from a biological point of view, 
it has been proved |33i that the exponential growth constant is strictly less than that of all self- 
avoiding loops, while the exponent (if it exists) is expected to coincide with that of all self-avoiding 
loops. 
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Fisher concluded that the above values of the loop class exponent c imply a continuous phase 
transition in c? = 2 and d = 3. We have to check the phase transition condition: In c? = 2, a 
numerical analysis of the oriented rooted SAP series data 19 gives xjj — 0.1436806285(8) and 
U{x^) = 0.6523866(2). The value of the generating function for straight paths is 1/(1 — xjj) < 
1.678. In d = 3, a corresponding analysis with unknotted oriented rooted SAP data and SAW data 
|28l App. C] gives xij = 0.045578(3) and U{x^) — 0.10(1). The path generating function value 
is bounded by that of random walks 1/(1 — Qxu) < 1.38. Thus, the phase transition condition is 
satisfied in both d = 2 and in d = 3 for straight paths. Note that in d = 2, self-avoiding walks 
as paths will result in no phase transition, since V{xu) > 2 for such paths, as follows from SAW 
data analysis. 

4 Self- Avoiding Poland-Scheraga models 

As discussed above, the previous PS models are not self- avoiding. Given a particular class of 
walks, for example, SAW, a PS model with segment-avoiding chains may be defined as follows. 
As paths, take only those walks with extremal first and last vertex: If v{0), f (1), . . . , v{n) are the 
vertices of an n-step walk v, this walk is taken as a path iff ^^^(O) < Vr^{i) < v^in) for all 1 < i < n. 
Such walks are bridges 28, Sec. 1.2], whose last step is in x-direction. For (unmarked) loops, 
take loops derived from the walks, as explained in Section im In our example these will be SAW 
loops. Marking of the loops may be achieved in different ways. To this end, consider for a given 
loop the sets Ai {Ar) of vertices of smallest (largest) x-coordinate. We distinguish four different 
types of marking: complete marking (with DNA constraint), where we mark a loop at all vertex 
pairs from Ai and Ar (whenever the DNA condition is satisfied), and unique marking (with DNA 
constraint), where we only mark a loop at a single vertex pair, for example the bottom vertex and 
the top vertex in a lexicographic ordering (if they satisfy the DNA condition). 

If the walk class is self-avoiding walks, this will result in self-avoiding PS models. Unique 
marking would then imply pn — 2p„ (we distinguish the two strands of a marked loop) and 
hence increase the previous exponents by one. Hence, such a PS model displays a first order 
transition with c = 5/2 in d = 2 and with c = 2.7631(18) in d = 3. (The phase transition 
condition is satisfied, as follows from the values U{x^) given in Section f3.3l and from the estimate 
V{xu) < 1/(1 — dxu) — 1 for bridges with last step in x-direction.) Unique marking with DNA 
constraint results in p„ < 2p„. If we assume that the exponential growth constant for marked 
self-avoiding loops is still given by fid (which we expect to be true, compare Section l5.2(l . this 
implies a critical exponent c greater or equal to the model with unique marking, i.e., a first order 
phase transition in d = 2 and d = 3. For complete marking, we have 2p„ < p„ < 2rtp„, which rules 
out a decrease of c by more than one. We however expect that the number of possible markings 
is of order 1 as n — + oo, such that c remains unchanged, and the model again displays a first 
order phase transition. Similar considerations apply for the case of complete marking with DNA 
constraint. 

5 Directed Poland-Scheraga models 

Wc present two classes of directed PS models which are exactly solvable and, by definition, take 
into account excluded volume interactions between all parts of the structure. The first class is 
derived from fully directed walks and is solvable in arbitrary dimension I15| . We will give explicit 
expressions for the generating function in d = 2 and in d = 3. This extends |25l I24ll^l27j . where 
the critical behavior has been derived. The model displays a first order phase transition only 
for d > 5. The second class, derived from partially directed walks, is considered in d = 2 only. 
Different variants of the second model display both a first order phase transition and a continuous 
phase transition in d = 2. 
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5.1 Fully directed walks and loops 

This directed PS model consists of fully directed walks for the paths in the chain. These only take 
steps in positive directions. The corresponding loops are staircase polygons, which consist of two 
fully directed walks, which do not intersect or touch, but have a common starting point and end 
point. Paths are attached to these points. We distinguish the two strands of a loop. This model 
satisfies all the assumptions discussed in Section [2.11 and also satisfies the DNA condition that 
the two segments of the loops are the same length. 

In (i = 2, the generating functions for paths and marked loops are |18) 

1 



V(x)^ , U(x) ^l-2x-Vl-ix. (5.1) 

1 ~ 2x 

U{x) is twice the generating function of staircase polygons. The loop class exponent is c = 3/2. 
We have xv — 1/2, xu — 1/4 and Wc = 1, i.e., a phase transition at T = oo. (If the empty path 
would not be allowed, a phase transition occurred at a finite temperature.) The free energy f{w) 
is given by 

/(i«)=log(-=i^=i^— 1 (l<u;<oo). (5.2) 
The fraction of shared bonds follows as 



n, ^ 2w; ^1 + (w; - 1)2 - 1 - (u; - 1)2/2 ,^ ^ , , , 

= 7-^ 7^ --^ ,^ , l<w;<oo, 5.3 

W-l 1 + (ly - 1)2 - ^1 + (w - 1)2 

which approaches zero linearly in w — 1. The asymptotic behavior of Zn{w) about w = 1 is given 
by 

Zn{w) ^ —^h^^/g^{w-l) (n ^ oo, w 1+), (5.4) 



n 



uniformly in w, where ha{x) is given by H2.15|l . 

The PS model of fully directed walks and loops is exactly solvable in arbitrary dimension jl5| . 
Consider fully directed walks on Z"*. If there are ki steps in direction i, the number of distinct 
walks, starting from the origin, is given by the multinomial coeflacient (''fc^l' ^ 'fc'''') • This results 
in Vd{x) — 1/(1 — dx). The generating function of a pair of such walks is given by 

The generating function Zd{x) can be interpreted as a chain, where each link consists of a staircase 
polygon or a double bond. Thus, it is related to the generating function Ud{x), where staircase 
polygons are counted twice, by 

Zd{x) = ^ . (5.6) 

I - [dx + Ud(x)) 

The functions Zd{x) satisfy Fuchsian differential equations of order d — 1, from which the singular 
behavior of Zd{x), and hence of Ud{x) can be derived jl5j . In dimensions 2 < d < 5, Ud{x) can be 
expressed in terms of Hcun functions. We get in d = 3 

^3(^)=(^) [l-Qx)-\l-x)-^K{k+)K{k^), (5.7) 
where K{k) is the complete elliptic integral of the first kind, and 

' \,x (^•«) 



Xi = 



9a;2 - 10a; + 1 
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We have xij = 1/9. In arbitrary dimension d > 3, it has been shown that the models have 
a critical point xjj ~ with exponent c = (d — l)/2, with logarithmic corrections in d = 3. 
Noting that Ud{x^) = 1 — 1/d— 1/Zd{x^) and Vd{xu) = d/{d— 1), we conclude that the phase 
transition condition Ud{x'ij)Vd{xu) < 1 is satisfied iff Zd{x^) < oo. This is the case in d > 4. 
For the corresponding PS model, the phase transition is thus first order for d > 6, and we have a 
continuous phase transition in dimensions 2 < d < 5, being at finite temperature in d > 4 only. 



5.2 Partially directed walks and loops 

We consider a directed PS model, where paths are partially directed walks on the square lattice 
Z^. These walks are self-avoiding and not allowed to take steps in the negative x direction. The 
number of paths of length n can be obtained 18 by considering the ways a walk of length n can 
be obtained from a walk of length n — 1. This yields the path generating function 

whose dominant singularity is a simple pole at xy = 1 — V^- The associated loop class is column- 
convex polygons, which has been analyzed in |7j, see Fig. ^ for an example. We assume unique 



Figure 1: A cohimn-convex polygon 



marking as explained in Section^ The (unmarked) loop generating function U{x) satisfies the 
algebraic equation of order four Tl, Eqn. 5.2] 

= - (a; - 1) V + (x^ - 7x^ + 3x - l){x - if U{x) 

+ 2(2x^ - Ux^ + lOx - 4)(x - i f U(xf 

' (5.10) 
-I- (5x3 _ 35^2 + 47^; _ 21) (x - 1) U{xf 

+ (2x3 _ 23^,2 33^ _ -^g^ U{xf. 

The dominant singularity at xjj = Xy = 3 — 2\/2 is a square-root singularity with an exponent 
c = 3/2. The critical amplitude is given by U{x^) = 0.086983(1). Thus 2J7(x^) V(x[/) < 1. This 
implies a continuous phase transition for the PS model at Wc = 1.922817(1). 

However, the loop class defined above does not obey the DNA constraint, which might increase 
the exponent c. Using standard techniques , it is possible to derive a functional equation for the 
generating function of column-convex polygons counted by perimeter and, in addition, by their 
upper and lower walk length. The corresponding expression is, however, much more difficult to 
handle than (|5.10|) , and there seems to be no practical way to extract the loop class exponent for 
the model with DNA constraint. On the other hand, it is possible to bound the loop class exponent 
from above by considering a refined (artificial) model. As we will argue below, this refined model 
has a loop class exponent c — 5/2. For the former model, we thus expect a value of 3/2 < c < 5/2, 
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with a value of c = 2 seeming plausible. This PS model would then be at the border between a 
first order and a continuous phase transition. 

For the refined model, let us further assume that the upper and lower walk in a loop are not 
allowed to touch the x-axis. The loop class of the refined model can be described in terms of 
bar-graph polygons |32| . These are column-convex polygons, one of whose walks is a straight line. 
Loops are composed of two bar-graph polygons, one of them in the half-plane y > 0, the other in 
the half-plane y < 0. Both polygons have the same horizontal coordinates and are constrained to 
have equal walk lengths. Paths of the model are partially directed walks, which are attached to 
the two loop vertices with vertical coordinate zero. 

Let (?m,n denote the number of bar-graph polygons with 2m horizontal and 2n vertical steps. 
The anisotropic perimeter generating function G{x,y) = gm,nX^^y"' satisfies an algebraic equa- 
tion of degree two j^., 

G{x, y}^xy+{y + + y))G{x, y) + xG^{x, y). (5.11) 

A generalization of the Lagrange inversion formula \\2' can be used to obtain a closed formula for 
its series coefficients. It is 

,) ^ - 1 ^ ("■: („ ,) („ ,) . ,..12, 

The number p2n of all combinations of bar-graph polygons subject to equal walk lengths is then 

(«-l)/2 

P2n = ^ 9l-2k,k- (5-13) 

The first few numbers P2n are 1,1,2,10,38,126,483,126,483,... for n = 3,4,5,.... We used 
standard methods of numerical series analysis |14j to estimate the critical point and critical expo- 
nent of U(x). An analysis with first order differential approximants, using the coefficients p2n for 
70 < n < 80, yields the estimates 

= 0.17157287(1), c = 2.5000(1) (5.14) 

This is, to numerical accuracy, indistinguishable from xj/ = 3 — 2\/2 ~ 0.171572875253..., as 
expected, and c = 5/2. Note that the phase transition condition 2U{x^)V{xu) < 1 is satisfied, 
since the loop generating function is bounded by that of column-convex polygons, which is known 
to satisfy the above inequality. 

6 Models extending the Poland- Scheraga class 

The PS model with oriented rooted self-avoiding loops reviewed in Section 1^131 displays a 
continuous phase transition in d = 2 and d = 2>. This led to the question whether introducing 
excluded volume effects between different segments of the chain can change the nature of the 
transition O EOl 121 ESI El El H • In Section^ we introduced self-avoiding PS models in order to 
answer that question in the affirmative, staying within the PS class. 

Another approach, which extends the PS class, consisted in simulations of self-avoiding walk 
pairs [HlElEini- We emphasize that these models are not PS models as defined in this paper, and 
that loop class exponents as defined in 12. 8|) bear no meaning for these models. The appropriate 
generalization is the exponent c' of the loop length distribution P{l,n) of the chain, which is 
assumed to behave in the hmit w ^ w+ like [201 EH E3 El El E| 

P{l,n) K r"' g{l/n) (1< ? < n, n cx)), (6.1) 

where g{x) is a scaling function, assumed to be constant in some analyses (201 El E3 El- Then, 
the same conclusions as in Section [2.21 aboiit the nature of the phase transition determined by the 
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loop length distribution exponent c' hold. This follows from the behavior of the mean loop length 
{l)n in the limit w — > w'^ , which is obtained from Ht).l|l by integration, 

n 

{l)^^^lP{l,n)f^ An^-""' (n^oo). (6.2) 

1=0 

If 1 < c' < 2, the mean loop length diverges, indicating a continuous transition. If c' > 2, the 
mean loop length stays finite, indicating a first order transition [SJH]. Note that assumption (|6.1|l 
is not fully consistent with PS models: The behavior 1)2. 23|) is consistent with (16.11) . which justifies 
the use of the same name for the two exponents. We found (|6.2(l to be satisfied for PS models if 
1 < c < 2 in Section ?IM If c > 2, however, the mean loop length exponent is independent of c 
for PS models — see (I2.24|) . This implies that a scaling form like (|6.1II cannot hold for PS models 
with loop class exponent c > 2. 

The simulation results arc c' = 2.44(6) ind — 2T and c' — 2.14(4) in d = 3 ^J, indicating a first 
order phase transition in d — 2 and in d = 3. This led to the conclusion that introducing excluded 
volume effects between different segments to the PS model of rooted self-avoiding loops drives the 
transition from continuous to first order |2(JI I21L I23| . The perturbative analytic arguments using 
the theory of polymer networks j20[ 1211 yield a very good approximation to the value of the 
loop length distribution exponent found in simulations. This suggests that the main mechanism 
responsible for the change of critical behavior in the PS model of oriented rooted self-avoiding 
loops is due to "local" excluded volume effects arising from forbidden attachments of paths to 
loops. 

7 Conclusion 

We discussed PS models with emphasis on the order of their phase transitions. We re-analyzed 
the old model of oriented rooted self-avoiding loops ^HI and found that the resulting PS model is 
not self-avoiding. Hence the conclusions about the order of the phase transition, which are most 
relevant to the recent discussions |21 120L 1211 1231 |S1 13 |22 , rely on a model where most 

chain configurations are not self-avoiding and thus cannot represent real DNA. Our self-avoiding 
PS model (unique marking, with self-avoiding bridges and (unrooted) self-avoiding loops as defined 
in Section 2J| yields a first order phase transition in both d ~ 2 and d = 3. It is expected that 
the other variants considered in Section ^ also yield a first order phase transition, but a detailed 
(numerical) analysis is needed to answer this question, as in the discussion in Section 0] 

For the presumably more realistic model of pairs of interacting self-avoiding walks |51 l2U[l?Tll23[ 
EliniQli our results suggest an interpretation of excluded volume effects, which complements the 
common one ^^2\^2^- The self- avoiding PS models of Section 21 correctly account for excluded 
volume effects within a loop, but overestimate excluded volume effects between different segments 
of a chain, due to their directed chain structure. Since this leads to a first order phase transition 
in d = 2 and c? = 3, we conclude that the relaxation of excluded volume effects between different 
segments of the chain does not change the nature of the transition. 

We also discussed several directed examples. One of these models displays a first order phase 
transition in c? = 2. Despite being exactly solvable, these models seem to be of limited relevance 
to the biological problem due to their directed structure and other limitations. 

The key question, as to which effects are responsible for the observed behavior in real DNA, 
is only partially answered by these results. With respect to directed models, which also account 
for different base pair sequences, it seems surprising however that simulations of melting curves 
|35[ |3j agree so well with experimental curves. In these simulations, a heuristic critical exponent 
accounting for the statistical weight of internal loops has to be inserted. It has recently been 
suggested that a value more realistic than the commonly used Fisher loop class exponent [TU| 
c ~ 1.763 might be the loop length distribution exponent of the polymer network approximation 
|2()j c' ~ 2.115. This suggestion implicitly assumes that the loop length distribution exponent of 
a non-directed model can be interpreted as an "effective" loop class exponent within a directed 
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model. Although appearing plausible, it seems difhcult to give this assumption a quantitative 
meaning. On the other hand, the MELTSIM simulation approach 4, is robust against a change 
of parameters. Indeed, melting curves can be simulated with satisfactory coincidence for both 
exponents, if the cooperativity parameter for the loops is adjusted accordingly 0. It would 
certainly be illuminating to obtain such simulated melting curves from a self-avoiding walk pair 
modeling ,2, 3 , 1 . We also mention the recent debate about the relevance of the self-avoiding 
walk pair model to real DNA [THII^ . 

The stiffness of the double stranded segments does not seem to qualitatively alter the critical 
properties of chains j8|, as can be seen by comparing the PS model result of Section 2.1, showing 
that critical properties are largely independent of the path generating function. The corresponding 
question of introducing energy costs for bending and torsion to single stranded segments has, to 
our knowledge, not been discussed. Even within the PS class, it is not clear how this affects the 
loop class exponents, since the changes just discussed might lead to loop classes with different 
exponential growth constants, but there is no argument as to how the exponents might change. 

In conclusion, the question of the mechanisms applying in real DNA which are responsible for 
the denaturation process and which explain multistep behavior as observed in melting curves, are 
still far from being satisfactorily answered in our opinion. 
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